#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>
#include <fstream>
#include "CubicPPSplines.h"

double f(double x) {
    return 1.0 / (1.0 + 25 * pow(x, 2));
}

int main() {
    double m_1, m_n;
    double epsilon = 1e-6;
    m_1 = (f(-1 + epsilon) - f(-1)) / epsilon;
    m_n = (f(1) - f(1 - epsilon)) / epsilon;

    std::vector<int> N = {6, 11, 21, 41, 81};
    std::vector<double> max_errors_complete, max_errors_natural, max_errors_periodic;

    for (int i = 0; i < N.size(); i++) {
        // Reset max errors for each N value
        double max_error_complete = 0.0;
        double max_error_natural = 0.0;
        double max_error_periodic = 0.0;

        std::vector<double> x, f_value;
        for (int j = 0; j < N[i]; j++) {
            x.push_back(-1.0 + 2.0 * j / (N[i] - 1));
            f_value.push_back(f(x.back()));
        }

        CubicPPSplinescomplete spline_complete(x, f_value, m_1, m_n);
        CubicPPSplinesnatural spline_natural(x, f_value, m_1, m_n);
        CubicPPSplinesperiodic spline_periodic(x, f_value);

        // Save data to file for plotting
        std::string dataFilename = "spline_data_N" + std::to_string(N[i]) + ".txt";
        std::ofstream outFile(dataFilename);
        int numDensePoints = 200;
        for (int k = 0; k < numDensePoints; k++) {
            double xi = -1.0 + 2.0 * k / (numDensePoints - 1);
            double y_true = f(xi);
            double y_complete_pp = spline_complete.value(xi);
            double y_natural_pp = spline_natural.value(xi);
            double y_periodic_pp = spline_periodic.value(xi);
            outFile << xi << " " << y_true << " " << y_complete_pp << " " << y_natural_pp << " " << y_periodic_pp << "\n";
        }
        outFile.close();

        // Calculate maximum error at midpoints of intervals
        for(int k = 0; k < N[i] - 1; k++){
            double x_middle_i = (x[k] + x[k + 1]) / 2.0;
            double error_complete = std::abs(spline_complete.value(x_middle_i) - f(x_middle_i));
            if (error_complete > max_error_complete) {
                max_error_complete = error_complete;
            }
            double error_natural = std::abs(spline_natural.value(x_middle_i) - f(x_middle_i));
            if (error_natural > max_error_natural) {
                max_error_natural = error_natural;
            }
            double error_periodic = std::abs(spline_periodic.value(x_middle_i) - f(x_middle_i));
            if (error_periodic > max_error_periodic) {
                max_error_periodic = error_periodic;
            }
        }

        // Store the maximum errors for this N
        max_errors_complete.push_back(max_error_complete);
        max_errors_natural.push_back(max_error_natural);
        max_errors_periodic.push_back(max_error_periodic);

        // Output the maximum errors for this N
        std::cout << "N=" << N[i] <<endl;
        std::cout<< " Maximum error for complete PP spline: " << max_error_complete << std::endl;
        std::cout << " Maximum error for natural PP spline: " << max_error_natural << std::endl;
        std::cout << " Maximum error for periodic PP spline: " << max_error_periodic << std::endl;

        // Generate gnuplot script and plot
        std::string scriptFilename = "plot_splines_N" + std::to_string(N[i]) + ".gp";
        std::ofstream scriptFile(scriptFilename);
        scriptFile << "set terminal pngcairo enhanced color size 1200,800\n";
        scriptFile << "set output 'spline_comparison_N" << N[i] << ".png'\n";
        scriptFile << "set title 'Comparison of Splines for N=" << N[i] << "'\n";
        scriptFile << "set xlabel 'x'\n";
        scriptFile << "set ylabel 'y'\n";
        scriptFile << "set key outside\n";
        scriptFile << "plot '" << dataFilename << "' using 1:2 with lines title 'True Function', \\\n";
        scriptFile << "     '" << dataFilename << "' using 1:3 with lines title 'Complete PP Spline', \\\n";
        scriptFile << "     '" << dataFilename << "' using 1:4 with lines title 'Natural PP Spline', \\\n";
        scriptFile << "     '" << dataFilename << "' using 1:5 with lines title 'Periodic PP Spline'\n";
        scriptFile.close();
        system(("gnuplot " + scriptFilename).c_str());
    }

    // Calculate convergence rates after all N values are processed
std::vector<double> convergence_rates_complete, convergence_rates_natural, convergence_rates_periodic;
for (size_t i = 1; i < N.size(); i++) {
    double rate_complete = -std::log(max_errors_complete[i - 1] / max_errors_complete[i]) / std::log(double(N[i - 1]) / double(N[i]));
    double rate_natural = -std::log(max_errors_natural[i - 1] / max_errors_natural[i]) / std::log(double(N[i - 1]) / double(N[i]));
    double rate_periodic = -std::log(max_errors_periodic[i - 1] / max_errors_periodic[i]) / std::log(double(N[i - 1]) / double(N[i]));
    convergence_rates_complete.push_back(rate_complete);
    convergence_rates_natural.push_back(rate_natural);
    convergence_rates_periodic.push_back(rate_periodic);
}

// Write convergence rates to file
std::ofstream convergenceFile("convergence_rates.txt");
for (size_t i = 1; i < N.size(); i++) { // Start from 1 because we have no previous value for N[0]
    convergenceFile << N[i] << " " << convergence_rates_complete[i-1] << " " << convergence_rates_natural[i-1] << " " << convergence_rates_periodic[i-1] << "\n";
}
convergenceFile.close();

// Plot convergence rates
std::string convergenceScript = "set terminal pngcairo enhanced color size 1200,800; "
                                "set output 'convergence_rates.png'; "
                                "set title 'Convergence Rates'; "
                                "set xlabel 'Number of Subintervals'; "
                                "set ylabel 'Convergence Rate'; "
                                "set key outside; "
                                "plot 'convergence_rates.txt' using 1:2 with lines title 'Complete', "
                                "'convergence_rates.txt' using 1:3 with lines title 'Natural', "
                                "'convergence_rates.txt' using 1:4 with lines title 'Periodic'";
system(("gnuplot -e \"" + convergenceScript + "\"").c_str());
    return 0;
}